%% Process function
function [TQ_plot,...
        TQ_tab,...
        cwt_struct] = ...
    process_fun(pstruct,...
        cur_trace,...
        TT_datetime_local,...
        batchin,...
        eventin)
    %% Calculation for T,Q, etc
    %josh crozier 2020
    try
        if ~isempty(cur_trace)
            %calculate T and Q
            [TQ_tab_local,cwt_struct,tqplot_temp] = ...
            T_Q_stack(cur_trace,...
                pstruct.Tmax,... %max period
                pstruct.Tmin,... %min period
                pstruct.minSTALTA,... %min amp/LTA ratio
                pstruct.minpromstd,... %min prominence of peaks (rel to std of data)
                pstruct.minampstd,... %min amp of peaks (rel to std of data)
                pstruct.minTseprat,... %minimum ratio of periods between peaks
                pstruct.mintsep,... %minimum separation between peaks in time
                pstruct.fitcycles,... %sets number of cycles to exp fit
                pstruct.minadjrsquare,... %min adjusted r squared value for exp fit
                pstruct.VoicesPerOctave,... %for cwt
                pstruct.WaveletWeights,... %for cwt: morse, amor, or bump
                pstruct.QWaveletWeights,...
                pstruct.WaveletParameters,... %morse wavelet params for cwt
                pstruct.NarrowWaveletParameters,...
                pstruct.specwindowtime,... %spectrogram window
                pstruct.decaythresh,... %min remaining amplitude fraction per cycle
                pstruct.Aweight,... %weight of mode "goodness" on amplitude
                pstruct.gweight,... %weight mode "goodness" on decay rate
                pstruct.rsquaredweight,... %weight mode "goodness" on exp fit r^2
                pstruct.promweight,... %weight mode "goodness" on prominence
                pstruct.startbounds,... %peak start time bounds
                pstruct.tqplot,... %plot
                pstruct.weight_decay_scale,... %decaying exp fit weights, smaller for faster decay
                pstruct.c_quantile,... %quantile at which to set constant offset for exp fits
                pstruct.c_std_frac,... %fraction of range beneath min to set c
                pstruct.ampcycles,... %number of cycles to sum over to check that is a local max
                pstruct.minhalflife,...%minumum oscillation half life
                pstruct.cwt_compress,... % output compressed cwt res {timewindow,numfreqbins} (empty for no compression)
                pstruct.fitstartoffsetcycles,...%offset to start expfit,recommend 1-2 cycles
                pstruct.LTAcycleoffset,... %amount to shift LTA (accouns for transform delay)
                pstruct.LTA_window,... %windowspanfor long-term avg
                pstruct.expmodelvars);  %variables to fit in "(d-c)*exp(bx)+c"

            if pstruct.tqplot
                TQ_plot = tqplot_temp;
            else
                TQ_plot = [];
            end

            if ~isempty(TQ_tab_local)
                %% calculate Q via FWHM method
%                 addpath('Fitspectra/Fitspectra/')
%                 data = cur_trace(1).data;
%                 save('fitspectradata.dat','data')
%                 fitspectra()
%                 addpath('QualityFactor')
%                 for jjj = 1:height(TQ_tab_local)
%                     startin = max([TQ_tab_local.startind(jjj)-...
%                         ceil(seconds(TQ_tab_local.T(jjj))*4*cur_trace(1).sampleRate), 1]);
%                     endin = min([TQ_tab_local.startind(jjj)+...
%                         ceil(seconds(TQ_tab_local.T(jjj))*20*cur_trace(1).sampleRate), cur_trace(1).sampleCount]);
%                     [pxx,f] = pspectrum(cur_trace(1).data(startin:endin),cur_trace(1).sampleRate);
%                     pxx = (pxx);
%                     for tracein_fwhm = 2:length(cur_trace)
%                         pxx = pxx + (pspectrum(cur_trace(tracein_fwhm).data,cur_trace(tracein_fwhm).sampleRate));
%                     end
%                     pxxtemp = pxx;
%                     Ttarget = seconds(TQ_tab_local.T(jjj));
%                     pxxtemp(1./f < Ttarget/2) = pxxtemp(1./f < Ttarget/2)/100;
%                     pxxtemp(1./f > Ttarget*4) = pxxtemp(1./f > Ttarget*4)/100;
%     %                 [f0, Q, params] = Zadler_et_al(f', pxxtemp', true)
%                     [f01, Q1] = full_width_half_max(f', pxxtemp');
% %                     [f02, Q2] = Robinson_and_Clegg(f', pxxtemp');
%                 end
                
                %% calculate Q via AR/Sompi method
%                 [Tsompi,Qsompi] = AR_Q_main(cur_trace,... traces:irisFetch trace struct
%                     TQ_tab_local.T,... %periods
%                     TQ_tab_local.start,...
%                     TQ_tab_local.start + 16*TQ_tab_local.T,...
%                     TQ_tab_local.Q);
%                 TQ_tab_local.Qsompi = Qsompi;
%                 TQ_tab_local.Tsompi = Tsompi;
                
                %% calculate amp and phases of ressonant periods
                ampphase_tab = table('Size',[height(TQ_tab_local),3*length(pstruct.channel_station_list)],...
                    'VariableTypes',repmat("double", 1, 3*length(pstruct.channel_station_list)),...
                    'VariableNames',["amp_" + pstruct.channel_station_list,...
                        "phase_" + pstruct.channel_station_list,...
                        "errphase_"+pstruct.channel_station_list]);
                ampphase_tab(:,:) = num2cell(NaN(height(ampphase_tab),width(ampphase_tab)));
                ampphase_tab = amp_phase(cur_trace, TQ_tab_local.T,...
                    TQ_tab_local.start,...
                    TQ_tab_local.start + pstruct.phasefitcycles*TQ_tab_local.T,...
                    pstruct.phase_WaveletParameters, ampphase_tab);
                %mean phase error
                ampphase_tab.mean_errphase = nanmean(ampphase_tab{:,"errphase_"+pstruct.channel_station_list},2);
                %mean phase error for best half of channels
                min50 = NaN(height(ampphase_tab),1);
                for kk = 1:height(ampphase_tab)
                    temp = ampphase_tab{kk,"errphase_"+pstruct.channel_station_list};
                    temp = temp(~isnan(temp));
                    temp = sort(temp);
                    temp = temp(1:ceil(length(temp)/2));
                    min50(kk) = mean(temp);
                end
                ampphase_tab.mean_min50_errphase = min50;
                TQ_tab_local = cat(2,TQ_tab_local,ampphase_tab);
                
%                 'plotting'
%                 plot_disp_map(TQ_tab_local(1,:))
                
                %% calculate first motions
                fm_tab = table('Size',[height(TQ_tab_local),length(pstruct.channel_station_list)+4],...
                    'VariableTypes',["datetime",repmat("double", 1, length(pstruct.channel_station_list)+3)],...
                    'VariableNames',["fm_start","fm_A_above_LTA","fm_STALTA","fm_Afrac",...
                        "fm_" + pstruct.channel_station_list]);
                fm_tab(:,1) = num2cell(NaT(height(fm_tab),1));
                fm_tab(:,2:end) = num2cell(NaN(height(fm_tab),width(fm_tab)-1));
                fm_tab = firstmotion(cur_trace, TQ_tab_local.T,...
                    TQ_tab_local.start - 6*TQ_tab_local.T,...
                    TQ_tab_local.start + 3*TQ_tab_local.T, ...
                    pstruct.minTseprat, pstruct.minSTALTA, pstruct.minampstd, ...
                    pstruct.fm_LTA_window_cycles,... %number cycles for STA/LTA first motion pick
                    pstruct.fm_minampfrac,... %min amp for first motion pick relative to max (prevents picking artifacts),...
                    fm_tab);
                TQ_tab_local = cat(2,TQ_tab_local, fm_tab);

                %calculate centroid

                %calculate attenuation

                TQ_tab = TQ_tab_local;
            else
                TQ_tab = table();
            end

            %save wavelet transform stack
            if pstruct.cwt_save
                if isempty(pstruct.cwt_compress)
                    cwt_savefun(strcat(pstruct.cwt_foldername,'/', ...
                        datestr(cur_trace(1).startTime, 30),'.mat'),...
                        TQ_tab_local, wt, time, s, coi)
                    cwt_struct = [];
                end
            end

        else
            TQ_plot = [];
            TQ_tab = table();
            cwt_struct = [];
        end %end calculations when trace not empty

    catch ME
        %check for TQ and etc errors
        disp(getReport(ME))
        disp(strcat('error TQ_stack:',num2str(batchin),',',...
            datestr(TT_datetime_local(eventin),31)))
        TQ_plot = [];
        TQ_tab = table();
        cwt_struct = [];
    end %end try

end %end function